#define INF     2e10f
#define DIMX     DIM_X
#define DIMY     DIM_Y
#define dx      (4.0f/DIMX)
#define dy      (2.0f/DIMY)
#define X0  -2.0f
#define Y0  -1.0f
#define Z0  -1.0f
#define ox 0.0f
#define oy 0.0f
#define oz 0.0f 
#define ENTITYS ENTITYS_NUM
#define DEPTH 50
#define ABS(x) (x>0.0f ? x:-x)
#ifndef MAX
#define MAX(a,b) (a > b ? a : b)
#endif
#ifndef MIN
#define MIN(a,b) (a > b ? b : a)
#endif
#ifndef M_PI
#define M_PI 3.1415f
#endif
// //////////////////产生随机数///////////////////////////////////////////////////
/*
Part of MWC64X by David Thomas, dt10@imperial.ac.uk
This is provided under BSD, full license is with the main package.
See http://www.doc.ic.ac.uk/~dt10/research
*/
#ifndef dt10_mwc64x_rng_cl
#define dt10_mwc64x_rng_cl

//! Represents the state of a particular generator
typedef struct{ uint x; uint c; } mwc64x_state_t;
typedef struct{ float3 O; float3 D; } ray;
enum{ MWC64X_A = 4294883355U };
enum{ MWC64X_M = 18446383549859758079UL };

void MWC64X_Step(mwc64x_state_t *s)
{
	uint X=s->x, C=s->c;
	
	uint Xn=MWC64X_A*X+C;
	uint carry=(uint)(Xn<C);				// The (Xn<C) will be zero or one for scalar
	uint Cn=mad_hi(MWC64X_A,X,carry);  
	
	s->x=Xn;
	s->c=Cn;
}
//! Return a 32-bit integer in the range [0..2^32)
uint MWC64X_NextUint(mwc64x_state_t *s)
{
	uint res=s->x ^ s->c;
	MWC64X_Step(s);
	return res;
}

#endif
inline float intersect_plane(float3 O, float3 D, float3 P,float3 N)
{
    // Return the distance from O to the intersection of the ray (O, D) with the
    // plane (P, N), or +inf if there is no intersection.
    // O and P are 3D points, D and N (normal) are normalized vectors.
    float denom = dot(D, N);
    if (ABS(denom) < 1e-6)
        return INF;
    float d = dot(P - O, N) / denom;
    if (d < 0)
        return INF;
    return d;
}

inline float intersect_sphere(float3 O,float3 D,float3 S,float R)
{
    // Return the distance from O to the intersection of the ray (O, D) with the
    // sphere (S, R), or +inf if there is no intersection.
    // O and S are 3D points, D (direction) is a normalized vector, R is a scalar.
	float3 SO=S-O;
    float b = dot(D, SO);
    float c = dot(SO, SO) - R * R;
    float disc = b * b - c;
    if (disc > 0 && b>0){
		float rt=sqrt(disc);
        return ((b-rt)>0 && (b+rt)>0)?(b-rt):(b+rt); //同时处理反射和折射
	}
    return INF;
}

inline float intersect(float3 O,float3 D,float4 obj,float3 normal,__global float4* debugarr)
{
	//半径无穷大，平面
    if (obj.w == INF)
	{
		// debugarr[0]=(float4)(INF,INF,INF,INF);
		//debugarr[0].xyz=normal;
		float t=intersect_plane(O, D, obj.xyz, normal);
		//debugarr[0].w=t;
		//debugarr[1].xyz=O;
		//debugarr[2].xyz=D;
        return intersect_plane(O, D, obj.xyz, normal);
	}
    else
        return intersect_sphere(O, D, obj.xyz, obj.w);
}

inline float find_intersect_object(float3 O,float3 D , int* obj_idx,__global const float4* pos,__global const float4* normal,__global float4* debugarr)
{
    float t = INF;
    *obj_idx = -1;
    for (int i=0;i< ENTITYS;i++){
		//debugarr[i]=normal[i];
		// debugarr[i].w=i+1;
        float t_obj = intersect(O, D, pos[i],normal[i].xyz,debugarr);
        if (t_obj < t){
            t= t_obj;
			*obj_idx =i;
		}
	}
    return t;
}

inline int find_shadow(float3 M,float3 N,float3 toL,int obj_idx,__global const float4* pos,__global const float4* normal,__global float4* debugarr)
{
	float ll=length(toL);
    for(int k=0;k< ENTITYS;k++){
        if (k != obj_idx){
            float t = intersect(M + N * 0.0001f, normalize(toL), pos[k],normal[k].xyz,debugarr);
            if(t<=ll) return 1;
		}
	}
    return 0;
}

inline float3 get_normal(int obj,float3 M ,__global const float4* pos,__global const float4* normal)
{
    // Find normal.
	float4 p=pos[obj];
    if (p.w==INF)
        return normal[obj].xyz;
    else{
		return normalize(M-p.xyz);
	}
}

inline float3 get_color(int obj,float3 M ,__global const float4* pos,__global const float4* color)
{
    // Find normal.
	float4 p=pos[obj];
    if (p.w==INF)
        return (float3)(((int)(M.x*2+2000)%2)==((int)(M.z*2+2000)%2) ? 0.9f : 0.1f);
    else{
		return color[obj].xyz;
	}
}

//计算反射光线
inline float3 reflect(float3 v, float3 n)
{
	return v - 2 * dot(v, n) * n;
}

// 计算折射光线
inline bool refract2(float3 v, float3 n, float ni_over_nt, float3* refracted)
{
	float3 uv = normalize(v);
	float dt = dot(uv, n);
	// 如果从法线负方向入射（例如球体内部的光线）
	if(dt>0){
		dt=-dt;
		n=-n;
		//ni_over_nt=1.0f/ni_over_nt;
	}
	else
		ni_over_nt=1.0f/ni_over_nt;
	float discriminant = 1.0 - ni_over_nt * ni_over_nt * (1 - dt * dt);
	if (discriminant > 0)
	{
		*refracted = ni_over_nt * (uv - n * dt) - n * sqrt(discriminant);
		return true;
	}
	else
		return false;
}

bool refract(const float3  v, const float3 n, float ni_over_nt, float3* refracted) {
    float3 uv =  normalize(v);
    float dt = dot(uv, n);
    float discriminant = 1.0 - ni_over_nt*ni_over_nt*(1-dt*dt);
    if (discriminant > 0) {
        *refracted = ni_over_nt*(uv - n*dt) - n*sqrt(discriminant);
        return true;
    }
    else
        return false;
}
inline float fresnel(float3 I, float3 N, float ior) 
{ 
    float cosi = dot(I,N); 
    float etai = 1.0f, etat = ior; 
    if (cosi > 0) 
	{ etai=ior;etat=1.0f ;} 
    // Compute sini using Snell's law
    float sint = etai/etat * sqrt(1.0f - cosi * cosi); 
    
    float cost = sqrt(1 - sint * sint); 
	cosi = ABS(cosi); 
	float Rs = ((etat * cosi) - (etai * cost)) / ((etat * cosi) + (etai * cost)); 
	float Rp = ((etai * cosi) - (etat * cost)) / ((etai * cosi) + (etat * cost)); 
	float kr = (Rs * Rs + Rp * Rp) / 2; 
    
    // As a consequence of the conservation of energy, transmittance is given by:
    // kt = 1 - kr;
	return kr;
} 
// 折射系数
inline float schlick(float cosine, float ref_idx) {
    float r0 = (1-ref_idx) / (1+ref_idx);
    r0 = r0*r0;
    return r0 + (1-r0)*pow((1 - cosine),5);
}

inline float3 random_in_unit_sphere(mwc64x_state_t *rng) {
    float3 p;
    do {
        p = 2.0f*(float3)((MWC64X_NextUint(rng)/4294967296.0f), (MWC64X_NextUint(rng)/4294967296.0f), (MWC64X_NextUint(rng)/4294967296.0f)) - (float3)(1.0f,1.0f,1.0f);
    } while (length(p) >= 1.0f);
    return p;
}

inline void savelight(float3 p1,float3 p2,int * counter,__global float4* debugarr,int i,int j)
{
	if(i==300 && j==137){
		debugarr[*counter*2].xyz=p1;
		debugarr[*counter*2+1].xyz=p2;
		*counter+=1;
		debugarr[0].w=*counter;
	}
}

// 产生0~1随机数
inline float RNG(mwc64x_state_t *rng)
{
	return (float)MWC64X_NextUint(rng)/4294967296.0f;
}

// 相机
float3 random_in_unit_disk(mwc64x_state_t *seed) {
    float3 p;
    do {
        p = 2.0f*((float3)(RNG(seed),RNG(seed),0.0f)) - (float3)(1.0f,1.0f,0.0f);
    } while (dot(p,p) >= 1.0f);
    return p;
}
inline ray camera_get_ray(float3 lookfrom, float3 lookat, float3 vup, float vfov, float aspect, float aperture, float focus_dist,float s,float t,mwc64x_state_t *seed)
 {
	// vfov is top to bottom in degrees
	float lens_radius = aperture / 2;
	float theta = vfov*M_PI/180;
	float half_height = tan(theta/2);
	float half_width = aspect * half_height;
	float3 origin = lookfrom;
	float3 w = normalize(lookfrom - lookat);
	float3 u = normalize(cross(vup, w));
	float3 v = cross(w, u);
	float3 lower_left_corner = origin  - half_width*focus_dist*u -half_height*focus_dist*v - focus_dist*w;
	float3 horizontal = 2*half_width*focus_dist*u;
	float3 vertical = 2*half_height*focus_dist*v;
	//ray get_ray(float s, float t)
	float3 rd = lens_radius*random_in_unit_disk(seed);
	float3 offset = u * rd.x + v * rd.y;
	origin=origin + offset;
	ray r={origin + offset, lower_left_corner + s*horizontal + t*vertical - origin - offset};
	return r;
}
__kernel void ray_tracing( __global const float4* color,__global const float4* pos,
						   __global const float4* normal,__global const float4* material,
						   __global const float4* lights, __global unsigned char * bitmap,
						   __global float4* debugarr
						   )
{
    // map from threadIdx/BlockIdx to pixel position
    int i = get_global_id(0);
    int j = get_global_id(1);
	if (i>=DIMX || j>=DIMY) return;
    int offset = i + j * get_global_size(0);
	// 随机采样 抗锯齿
	mwc64x_state_t seed={i*9000+19797,j*1040+1432};
	MWC64X_NextUint(&seed);
    float3 O =(float3)(ox,oy,oz); //眼睛的位置
	
	int counter=0;//光线计数
	int obj_idx=0;// 与视线相交的实体的编号

	float3 color_ij=(float3)(0,0,0);
	
	int ns=100;//采样次数
	//float3 tmp=(float3)(0.0f);
	for (int s=0;s<ns;s++){
		float reflection=1.0f;
		float3 tmp=(float3)(1.0f);
		// 相机
		float3 lookfrom=(float3)(0,0,0);
		float3 lookat=(float3)(0,0,-20);
		float dist_to_focus = length(lookfrom-lookat);
		float aperture = 0.0;
		// float3 lookfrom=(float3)(13,2,3);
		// float3 lookat=(float3)(0,0,0);
		// float dist_to_focus = 10.0;
		// float aperture = 0.0f;
		float u = (float)(i + RNG(&seed)) / (float)(DIMX);
        float v = (float)(j + RNG(&seed)) / (float)(DIMY);
		ray r=camera_get_ray(lookfrom,lookat, (float3)(0,1,0),30,(float)DIMX/(float)DIMY,aperture,dist_to_focus,u,v,&seed);
		//float3 Q =r.D; //视线和屏幕的交点
		//float3 Q =(float3)(X0+(i)*dx,Y0+(j)*dy,Z0); //视线和屏幕的交点
		float3 D=normalize(r.D);// 视线的方向
		float3 ORG=r.O;
		for (int n=0;n<DEPTH;n++){
			float dt= find_intersect_object(ORG, D,&obj_idx,pos,normal,debugarr);
			if( dt == INF) {// 不相交直接退出
			float t = 0.5f*(D.y + 1.0f);
			if (n==DEPTH-1)
				tmp*=0;// 光线被困住，全黑
			else
				tmp*=((1.0f-t)*(float3)(1.0f, 1.0f, 1.0f) + t*(float3)(0.5f, 0.7f, 1.0f));//背景色
			break;
			}else{
			// 交点
			float3 M = ORG + D * dt;
			// 保存视线
			//savelight(ORG,M,&counter,debugarr, i, j);
			// 法线
			float3 N = get_normal(obj_idx,M ,pos,normal);
			// // 指向光源
			// float3 L=lights[0].xyz;
			// float3 toL =L - M;
			// // 指向眼睛
			// float3 toO = normalize(O - M);
			// // 判断光源是否被挡住
			// int shadowed=find_shadow(M,N,toL,obj_idx,pos,normal,debugarr);
			// toL=normalize(toL);
			// // 材质参数
			float4 mat=material[obj_idx];
			// 颜色属性
			float3 mycolor=get_color(obj_idx, M ,pos,color);
			tmp*=mycolor; //或者其他方法混合颜色
			// // 计算环境光
			// temp+=0.05f;
			// if (!shadowed){
			// // 计算漫反射 Lambert shading (diffuse).
				// temp+=mat.y*MAX(dot(N, toL), 0) * get_color(obj_idx,M ,pos,color);
			// // 计算镜面反射 Blinn-Phong shading (specular).
				// temp+=mat.z*pow(MAX(dot(N, normalize(toL + toO)), 0), 50) * lights[1].xyz;
			// }
			// 按材质计算将要追踪的下一条反射光线
			if(mat.x==0.0f)
			{ // 纯漫反射表面
				D=normalize(N+random_in_unit_sphere(&seed));
			}
			else if(mat.x==1.0f)
			{ // 纯反射表面
				if (mat.w==INF)
				{  	// 不透明
					float fuzz =0.0f; //模糊因子
					D=normalize(reflect(D, N)+fuzz*random_in_unit_sphere(&seed));
				}
				else //透明，计算折射光线
				{	float3 outward_normal;
					float ref_idx,ni_over_nt,reflect_prob,cosine;
					ref_idx=mat.w;
					if (dot(D, N) > 0){//光线在玻璃内部
					outward_normal -=N;
					ni_over_nt = mat.w;
					
					// cosine = ref_idx * dot(r_in.direction(), rec.normal) / r_in.direction().length();
					cosine = dot(D, N);
					cosine = sqrt(1 - ref_idx*ref_idx*(1-cosine*cosine));
					}
					else {
						outward_normal =N;
						ni_over_nt = 1.0f / ref_idx;
						cosine = -dot(D, N);
					}
			 
					float3 refrac; //折射光线
					//if(refract(D,outward_normal,ni_over_nt,&refrac)) //可以折射
					if(refract2(D,N,ref_idx,&refrac)) {//可以折射
					 
					reflect_prob = schlick(cosine, ref_idx);//反射的概率}
					reflect_prob=fresnel(D,  N,  mat.w);
					}	// { 
							// if (RNG(&seed)<0.5) 
								//D=refrac;
							// else
								// D=normalize(reflect(D, N));
								
						// }
					else {
						//D=normalize(reflect(D, N));
						reflect_prob=1.0f;
					}
					if(RNG(&seed)<reflect_prob)
						 D=normalize(reflect(D, N));
					else
						 D=refrac;
				}
			}
			ORG=M+0.0001f*dot(N,D)*N; //反射光线的起点
			}
		}
		color_ij+=tmp;
	}
	color_ij/=ns;
	color_ij.x=sqrt(color_ij.x);
	color_ij.y=sqrt(color_ij.y);
	color_ij.z=sqrt(color_ij.z);
    bitmap[offset*4 + 0] = (int)(MIN(color_ij.x,1) * 255);
    bitmap[offset*4 + 1] = (int)(MIN(color_ij.y,1) * 255);
    bitmap[offset*4 + 2] = (int)(MIN(color_ij.z,1)  * 255);
    bitmap[offset*4 + 3] = 255;
}